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Abstract 

A new algorithm for the derivation of low-density series for percolation on directed 
lattices is introduced and applied to the square lattice bond and site problems. Nu- 
merical evidence shows that the computational complexity grows exponentially, but 
with a growth factor A < \/^, which is much smaller than the growth factor A = v^ of 
the previous best algorithm. For bond (site) percolation on the directed square lattice 
the series has been extended to order 171 (158). Analysis of the series yields sharper 
estimates of the critical points and exponents. 

1 Introduction 

Directed percolation (DP) [|l], ^ can be thought of as a purely geometric model in which the 
bonds or sites of a hyper-cubic lattice Z'^ are either present with probability p or absent (with 
probability q = 1 — p). Unless otherwise specified I shall be looking at bond percolation on 
the directed square lattice. As in ordinary percolation one might be interested in the average 
cluster-size S{p). The difference is that in DP connections are only allowed along a preferred 
direction given by an orientation of the edges. Fig. ^ shows the part of the square lattice 
which can be reached from the origin using no more than five steps. Apart from the inherent 
theoretical interest DP has been associated with a wide variety of physical phenomena. In 
static interpretations, the preferred direction is spatial, and DP could model gravity driven 
percolation of fluid through porous rock with a certain fraction of the channels blocked 
0, crack propagation Q or electric current in a diluted diode network 0. In dynamical 
interpretations, the preferred direction is time and one realization is as a model of an epidemic 
without immunization. DP type transitions are also encountered in many other situations 
including Reggeon field theory [§], chemical reactions 0, heterogeneous catalysis and other 
surface reactions [Q , self-organized criticality , and even galactic evolution []TU| . 



Domany and Kinzel [TT|] demonstrated that bond percolation on the directed square 
lattice is a special case of a one- dimensional stochastic cellular automaton. The evolution 
of the model is governed by the transition probabilities W{ax\(Ji, cfr), with o"i = 1 if site i is 
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occupied and otherwise. It is the probabihty of the site x being in state ax at time t given 
that the sites x — 1 and x + 1 at time t — 1 were in states o"; and 0"^, respectively. One has 
a free hand in choosing the transition probabihties as long as one respects conservation of 
probability, P{l\ai,ar) = 1 — P(0|cr;, cr^). Bond percolation corresponds to the choice: 

W{0\ai,ar) = {l-pr^''^ (1) 

Note that while one talks of occupied sites in the cellular automata language, the choice of 
W above does indeed correspond to bond percolation in the static version, as can easily be 
verified by explicitly looking at all possible configurations and their associated probabilities. 
The behavior of the model is controlled by the branching probability (or density of bonds) 
p. When p is smaller than a critical value Pc all clusters remain finite. Above Pc there is a 
non-zero probability of finding an infinite cluster and S{p) diverges as p -^ p~ . In the low- 
density phase {p < Pc) many quantities of interest can be derived from the pair-connectedness 
Cx,t{p), which is the probability that the site x is occupied at time t given that the origin 
was occupied at t = 0. The moments of the pair-connectedness may be written as 

oo 

/^n,n^(p)=5Z5^x"ra,t(p). (2) 

t=0 X 

In particular the average cluster size S{p) = /io,o(p)- Due to symmetry, moments involving 
odd powers of x vanish. The remaining moments diverge as p approaches the critical point 
from below: 

^^in,miP) OC (Pc - pY^^+^-^^^-\\\ P^P~. (3) 

For isotropic percolation many exact results are known for two-dimensional systems, e.g., 
the critical exponents are known exactly from arguments based on conformal field theory 
]13| and for some lattices the critical point is known exactly |]12| . Isotropic bond percolation 



is also related to the g ^ 1 limit of the Potts model |]14|, [T^. For directed percolation no 
such exact results are known, except for special limited versions such as compact directed 



percolation [16]. Though it is worth mentioning that directed bond percolation is related to 



the g — !► 1 limit of the chiral Potts model |T^ and to the m — > limit of the m-friendly walker 



problem |jT8|, 0. So in order to study directed percolation one has to resort to numerical 
methods. For many problems the method of series expansions is by far the most powerful 
method of approximation. For other problems Monte Carlo methods are superior. For the 
study of DP on two-dimensional lattices, series analysis is undoubtedly the most appropriate 
choice. This method consists of calculating the first coefficients in the expansion of fin,m{p)- 
Given such a series, using the numerical technique known as differential approximants ||20|| , 
highly accurate estimates can frequently be obtained for the critical point and exponents, as 
well as the location and critical exponents of possible non-physical singularities. 

The difficulty in the enumeration of most interesting lattice problems is that, compu- 
tationally, they are of exponential complexity. Initial efforts at computer enumeration of 
square lattice directed percolation were based on direct counting. The computational com- 
plexity is proportional to A", where n is the number of terms in the series, and Ai = 3, is the 



connective constant for directed lattice animals on the square lattice [^. A dramatic im- 
provement was achieved by the transfer matrix technique of Blease who devised an algorithm 
with a complexity, which is proportional to A2, where A2 = a/2 ca 1.414. This work resulted 
in a derivation of the first 21 terms in the low-density series on the square lattice. Over the 
years further significant progress was obtained by Essam and co-workers who supplemented 



the transfer matrix technique by a weak subgraph expansion and obtained first 25 terms |23 



and some years later 35 terms |2^. Note that while the weak subgraph expansion enables 
one to derive more terms it does not improve the computational complexity. Eventually they 
devised a resummation technique |2^ which was used to derive the so-called non-nodal graph 



expansion [^. This in turn enabled them to obtain twice as many terms [^ as could be 
obtained by the bare algorithm of Blease, resulting in a series of 49 terms. So this approach 
has a complexity, which is proportional to Ag, where A3 = \f2 ~ 1.189. This method was 



later supplemented by an extension method |23] , originally suggested by work of Baxter and 



Guttmann |^, based on predicting correction terms from successive calculations on finite 
lattices of increasing size. Again this extension method does not improve the computa- 
tional complexity but it did, in conjunction with several years of improvements to computer 
technology, result in an extension of the series to 112 terms. In this paper I shall describe 
a new algorithm based on a direct calculation of the non-nodal graph expansion counting 
only graphs with non-zero weights. This has reduced both time and storage requirements 
by virtue of a complexity proportional to A4, where the numerical evidence suggest that 
A4 < \pl ~ 1.091. The series has now been calculated to 171 terms using essentially the 
same computational resources used earlier to derive 112 terms, but without the need for a 
cumbersome and complicated extension procedure. The series for the directed square lattice 
site problem has been extended to 158 terms from the previous best of 106 terms. 

In the next section I will briefiy review the previous best method for deriving low-density 
series expansion for square lattice directed percolation. In section ^ I give a detailed de- 
scription and empirical analysis of the computational complexity of the improved algorithm. 
The results of the analysis of the series are presented in Section ^. 

2 Series expansions 

From Eq. (^ it follows that the first and second moments can be derived from the quantities 

S{t) = Y,C.Av) and X(t) = J]x2a,t(p). (4) 

X X 

S{t) and X{t) are polynomials in p obtained by summing the pair-connectedness over all 
lattice sites whose parallel distance from the origin is t. It has been shown ^9] that the 
pair-connectedness can be expressed as a sum over all graphs (or finite clusters) formed by 
taking unions of directed paths connecting the origin to the site {x,t), 

a,b) = $^rf(<7)pi^i, (5) 

g 

where l^^l is the number of bonds in the graph g. Any directed path to a site whose parallel 
distance from the origin is t contains at least t bonds. From this it follows that if S{t) and 
X{t) have been calculated for t < t^g,^ then one can determine the moments to order tmax- 
One can however do much better as demonstrated by Essam et al. [E3]. They used a non- 



nodal graph expansion, based on work by Bhatti and Essam ||2^ , to extend the series to order 
2tmax + 1- The non- nodal graph expansion has been described in detail in |2^ and here I will 
only summarize the main points and introduce some notation. A graph g is nodal if there is 
a point (other than the terminal point) through which all paths pass. It is clear that each 
such nodal point effectively works as a new origin for the cluster growth. This is the essential 
idea behind the non- nodal graph expansion. Let S^{t) and X^{t) be the contribution from 
non-nodal graphs to S{t) and X{t), respectively. The non- nodal expansions are obtained 



recursively from the polynomials S{t) and X{t) [|25]. Next one forms the moments S^ , fi^i, 
11q2) ^^d /i^Q of non- nodal contributions equivalent to Eq. (0). The final series are obtained 
from the formulas 

(6) 

(7) 



S -- 


= 1/(1 -S^) 


1^0,1 = 


= <1^^ 


/^0,2 = 


= [<2 + 2Ki)^^]^^ 


1^2,0 = 


- 11^ 92 



(9) 

Further extensions of the series can be obtained by using a procedure similar to that of 
Baxter and Guttmann p^. One looks at correction terms to the series and tries to identify 



extrapolation formulas for the first n^ correction terms allowing one to derive a further n^ 



series terms correctly. Details of this procedure can be found in p7 |. 



The pair-connectedness can be calculated ^^ using a transfer-matrix technique based 



on the Domany-Kinzel model. The probability of finding a given configuration of sites at 
time t was calculated by moving a boundary through the lattice one site at a time. Fig. ^ 
shows how the boundary (marked by large filled circles) is moved in order to pick up the 
weight associated with a given 'face' of the lattice at a position x along the boundary line. 
At any given stage this hne cuts through a number of, say k, lattice sites thus leading to a 
total of 2^ possible configurations along this line. Configurations along the boundary line are 
trivially represented as binary numbers, and the probability of each configuration is given by 
a truncated polynomial in p. Let 5*0 = (ai, ...., 0"^_i, 0, cr^+i, ..., 0"^) be the configuration of 
sites along the boundary with at position x and similarly 5*1 = (cti, ...., (Xx^i, 1, cTx+i, ■■■, crk) 
the configuration with 1 at position x. Then in moving the boundary at the a;'th site the 
polynomials are updated as follows 

P{SO) = W{0\0,ai)PiSO) + W{0\l,ai)PiSl), 
P{S1) = W{l\0,ai)P{SO) + W{l\l,ai)P{Sl). 

The pair-connectedness is obviously symmetrical in x, Cx,t{p) = C-x,t{p), so it suffices 
to calculate the pair-connectedness for x > 0. More importantly, due to the directedness of 
the lattices, if one looks at sites (x, t) with x > they can never be reached from points 
{x', t') in the part of the lattice for which t' > [t/2\ , x' < — [t/2j . In Fig. |l| the part of the 
lattice needed in order to calculate C{x,t) up to t = 5 is enclosed by the dashed lines and 
the two sites on the boundary-line to the right of the dashed line can be disregarded. In 
more general terms this means that the pair-connectedness at points a distance t from the 
origin can be calculated using a boundary which cuts through at most [t/2j + 1 sites. Thus 
the memory (and time) required to calculate S{t) and X(t) grows like 2l-*/^J+^. Using the 
non-nodal expansion it follows that the time and memory required to calculate the series to 
a given order n grow like 2"/^, and thus that the computational complexity of this algorithm 
is exponential with growth-factor A = v2- In |^ this calculation was performed up to 



^max = 47. By looking at correction terms the series for S{p), Ho,i{p) and /io,2(p) were 
extended to order 112 and the series for fi2,o{p) to order 111. 

3 The new algorithm 

In the following I shall describe a new algorithm based on calculating the pair-connectedness 
using the expression in terms of contributing graphs given in Eq. (^ . The algorithm directly 



calculates the non-nodal contributions and the numerical evidence shows that asymptotically 
the computational complexity has a growth-factor A < \/2. As already noted, the directed 
lattice weights d{g) are non-zero if and only if the graph g is a union of paths from the origin 
to {x,t). The graphs with non-zero weights form a very limited set of all possible graphs 
passing through {x,t). The restriction to unions of paths is very strong and one immediate 
consequence is that graphs with dangling/dead-end parts make no contribution. Likewise a 
graph without dangling parts makes a contribution to C{x,t) only if it terminates exactly 
at (x, t). It could of course contribute at a later stage when the various branches join. Fig. § 
illustrates these points. The thick solid lines shows an example of a graph contributing to 
C{x, t). The dotted lines illustrate dangling ends which are not allowed, and the thin dashed 
lines show branches which would have to be joined at a later stage in order for this extended 
graph to make a contribution to C{x',t'). 

The directed weights have the further properties [^ 

d{g)=dig')d{g"), rf(^) = 0, ±1. (10) 

where g' and g" are subgraphs of g. So the d-weights factorize and for contributing graphs 
are either 1 or —1. Furthermore, for contributing graphs d{g) = (— 1)*^^)+^, where t{g) is 
the number of independent paths from the origin to the end point of the graph g. The 
calculation of d{g) is in principle quite complicated were it not for the factorization property 
and the fact that factors of —1 are picked up each time two paths join (or alternatively when 
one path splits into two paths). Indeed one sees that t{g) is simply the number of times 
a path splits into two, which of course equals the number of times two paths join to form 
a single path. This means that once again the pair-connectedness can be calculated via a 
transfer-matrix type algorithm by moving a boundary line through the lattice one row at 
a time with each row built up one site at a time. The sum over all contributing graphs is 
calculated as one goes along. 

In Fig. I have given a pictorial representation of how the boundary line polynomials are 
updated by the new algorithm. First it should be noted that there are two states per site 
with the following prescription: cxj = 1 if a bond has been inserted along an edge from the 
row above, and o"j = otherwise. We are looking at a move similar to that in Fig. |I|, though 
in this case the sites involved in the move are the one on the top and its neighbours in the 
row below. I shall refer to the configuration after the move as the 'target' configuration and 
it accounts for the state of the sites on the bottom row. The configurations prior to the 
move, given by the state of the top and bottom right sites, shall be called 'source' states. 
The state of sites away from the kink in the boundary does not influence the updating. In 
Fig. 1^ sites with incoming bonds are marked with a filled circle, those with open circles have 
no incoming bonds, sites marked by a shaded circle have incoming bonds in the target state 
but not in the source state, and bonds are marked by thick lines. Note that the avoidance 
of dangling ends is easily achieved by ensuring that sites with incoming bonds have at least 
one outgoing bond. In the first equation the target configuration, a,,, has bonds entering 
both sites and is fed by the source configurations a.^o and a,^,. In moving the boundary 
from the top site, bonds can be inserted to the left and/or right. The source state cr.^o has 
no bonds entering the left site so both bonds have to be added with an associated weight p^. 
From the source state a,^, a bond has to be inserted on the left edge, but on the right edge 
a bond can be either absent, in which case an associated weight p is required, or a bond can 
be added giving a weight —p^ because two bonds were inserted and two paths join on the 
site to the right. Similar considerations lead to the other equations. 

Limiting the calculation to non-nodal contributions is very simple; whenever the bound- 
ary line reaches the horizontal position all one need do is set to zero the polynomials of states 



with a single incoming bond, i.e., states represented by integers of the form 2^ . This obvi- 
ously ensures that configurations with a nodal point are deleted from the calculation. The 
pair-connectedness at the following time can be calculated from the states with incoming 
bonds at nearest neighbour sites and no incoming bonds on any other sites. 

In a calculation to a given order n we need to calculate the non-nodal contributions 
for all t < tinax = n/2. For a given t' < t the possible configurations along the boundary 
line is limited by constraints arising from the facts that graphs have to terminate at level 
t and have no dangling parts. As mentioned, the "no dangling parts" restraint is equiva- 
lent to demanding that sites with incoming bonds also have outgoing bonds. Therefore a 
configuration for which the maximal separation between sites with incoming bonds is r will 
take at least another r — 1 steps before collapsing to a configuration with a pair of nearest 
neighbour sites with incoming bonds. Consequently if t' + r > t that configuration makes no 
contribution to C(x, t) for any x and can be discarded. Furthermore the minimum number 
of bonds needed to build up the internal structure of the configuration (that is the occupied 
sites except for the left- and right-most sites) will have to be inserted during the collapse of 
the configuration. Since we are interested only in non-nodal graphs we further know that we 
have to insert at least 2t bonds to get from the origin to a point (x, t). It is easy to calculate 
the minimum order, A^min, of the boundary polynomial as the configuration is built up, but 
in so doing we have also counted 2t' bonds coming from the paths leading to the left- and 
right-most sites, respectively. So the minimal order to which a configuration contributes, 
Ncont, from row t' is 

Accent = 2A^^in + 2t - 4t'. (11) 

If Acont > 1^ the configuration can be discarded since it will only contribute at an order 
exceeding that to which we want to carry out our calculation. Further memory savings are 
obtained by observing that in calculating C{x, t) we know that the non-nodal graphs have 
at least 2t bonds, so we need only store n — 2t coefficients, and when the boundary is moved 
from one row to the next we discard the two lowest order terms in the boundary polynomials 
since they are all zero. 

The algorithm for calculating the series for directed site percolation is very similar, but of 
course we have different rules for updating the boundary polynomials. Though it should be 
noted that the meaning of 'O's and 'I's in the encoding of configurations is changed slightly 
to signify unoccupied and occupied sites respectively. The updating rules are easy to derive 
from those for the bond problem. In fact all we need to do is take the rules depicted in 
Fig. ^ and wherever we have inserted a bond on the right to a site already occupied the 
weight associated with that possibility is simply divided by p because in site percolation we 
cannot reoccupy a site already occupied. This leads to the following update rules for the 
site problem: 

P(a.,.) = p'P(a.,o), 

P(ao,.) = pP(a.,o)+P(ao,.)-P(a.,.), 

P(a.,o) = pP(a.,o), (12) 

P(ao,o) = P(ao,o). 

The only other change is that the formula for the minimal order of a contributing configu- 
ration is changed a little, becoming 

Ar,o„t = 2N^i^ + 2t- Nocc - ^t' (13) 



where we have subtracted the number of occupied sites, A^occ, in the configuration (otherwise 
they would have been counted twice from the first factor). 

In Fig. ^ I have plotted the number of distinct configurations which make a contribution 
in a calculation of the square lattice bond series to order n. As can be seen the numeri- 
cal evidence clearly shows that the asymptotic growth in the number of configurations is 
exponential with a growth factor A < \/2. As already noted, this is a very substantial 
improvement on the previous best algorithm which had a growth factor A = \/2. In plain 
words this means that whenever computer memory is doubled one can derive an additional 
8 terms with the new algorithm as compared to an additional 4 terms with the previous 
best algorithm. Using all the memory minimization tricks mentioned above it is possible 
to derive the series for moments of the pair-connectedness to order 171 with a maximum of 
approximately 5.5Gb of memory. The old algorithm would have required more than 250 000 
times as much memory and even if the extension procedure was used to derive say the last 
20 terms correctly more than 10 000 times as much memory would have been required. For 
directed site percolation the series has been derived to order 158 using similar computational 
resources. 

Finally a few remarks of a more technical nature. The number of contributing configura- 
tions becomes very sparse in the total set of possible states along the boundary line and as 



is standard in such cases one uses a hash-addressing scheme [^. It is probably also worth 
mentioning that since the update involves only two nearest neighbour sites and does not 
depend on the state of other sites along the boundary, this algorithm lends itself very nat- 
urally to parallel computing. Since the integer coefficients occurring in the series expansion 



become very large, the calculation was performed using modular arithmetic [32|. This in- 
volves performing the calculation modulo various prime numbers pi and then reconstructing 
the full integer coefficients at the end. In order to save memory I used primes of the form 
Pi = 2^^ — Tj so that the residues of the coefficients in the polynomials can be stored using 
16-bit integers. The Chinese remainder theorem ensures that any integer has a unique repre- 
sentation in terms of residues. If the largest absolute values occurring in the final expansion 
is m, then we have to use a number of primes k such that piP2 ■ ■ ■Pk/'^ > ^^- Note that it 
is not necessary to be able to uniquely reproduce the intermediate values, which in some 
cases can be much larger than the final ones. Up to 12 primes were needed to represent the 
coefficients correctly. 

The updating rules given above are quite general and can be used for calculating the pair- 
connectedness of percolation problems on the directed square lattice in many other cases. 



One example is that of a lattice with an impenetrable wall [^, as shall be demonstrated in a 
forthcoming article, though in this case the calculation is not confined to non-nodal graphs. 
Other applications include directed percolation with temporal or spatial disorder in which 
the branching probability p depends on t or a; |^ . The basis for the algorithm used in this 
paper is Eq. (^ . It is a very general expression for the pair connectedness which in fact should 
hold on any directed lattice and also for more complicated processes. Thus the work started 
in this paper can be generalized to other planar lattices or to higher dimensional lattices, 
though naturally the updating rules and other details of the actual implementation will 
vary from application to application. Another very interesting possibility is the application 
to unidirectionally coupled directed percolation [^ ^, |3^. In this case one studies a 



hierarchy of identical directed percolation processes. The O'th level process is just an ordinary 
directed percolation process. The fc'th level process is a directed percolation process evolving 
according to the usual rules. But in addition sites (x, t) may become occupied at some given 
rate if the corresponding site was occupied at level {k — 1). Studies of this system showed that 
the exponents depend on /c, in essence showing that as k is increased percolation becomes 



easier. By truncating the hierarchy at some low value of k it should be possible to generalize 
the algorithm to study this new and very interesting problem. 

4 Analysis of series 

The series for moments of the pair-connectedness were analyzed using differential approx- 
imants. A comprehensive review of these and other techniques for series analysis may be 
found in |2^. This allows one to locate the critical point and estimate the associated critical 
exponents fairly accurately, even in cases where there are additional non-physical singular- 
ities. Here it suffices to say that a i^th-order differential approximant to a function /, for 
which one has derived a series expansion, is formed by matching the coefficients in the poly- 
nomials Qi and P of order iVj and L, respective, so that the solution to the inhomogeneous 
differential equation 

j2Q^i^)i^^yfi^) = pi^) (14) 

agrees with the first series coefficients of /. The equations are readily solved as long as the 
total number of unknown coefficients in the polynomials is smaller than the order of the 
series n. The possible singularities of the series appear as the zeros Xi of the polynomial Qk 
and the associated critical exponent Aj is estimated from the indicial equation 



The physical critical point is generally the first singularity on the positive real axis. 

4.1 The bond series 

In this section I will give a detailed account of the results of the analysis of the square 
bond series. The analysis of the square site is described in the following section. In addi- 
tion to the moment series I have also analyzed the series /io,2(p)/Aio,i(p) ~ (Pc — p)'"^^ and 
/i2,o(p)/^o,2(p)/(/xo,i(p))^ ~ (Pc-p)'^"^. 

As one increases the order Ni and L of the polynomials in Eq. |14| and thus the number of 
terms used to form the differential approximants, one would generally expect to obtain more 
accurate estimates for the critical parameters. Likewise one often finds that the estimates 
show some dependence on the order K of the differential approximant and/or the order 
L of the inhomogeneous polynomial. In previous work [^ it was observed that for some 



series the estimates from first-order differential approximants showed a marked change with 
increasing L. Analysis of the longer series of course confirm this and also that some series 
(in particular /io,i(p) and /io,2(p)) show a certain drift as K is increased. But over all the 
estimates from the various series are exceptionally well converged and for K >2 show little 
if any change as L is changed or K increased. 

In order to gauge the effect of any systematic drift and lack of convergence to the true 
critical values it is helpful to plot the estimates of Pc and associated exponents obtained from 
approximants to the various series as a function of the number of terms. Plotted in Fig. ^ 
are estimates for pc obtained from third-order differential approximants with L increasing 
from to 50 in steps of 5. The orders Ni were chosen so that the difference between 
the order of the polynomials Qi never exceeded 1. Each point in the figure represents an 
estimate from a single approximant. From this study, two conclusions are immediately 



obvious. First of all most of the series show some drift in the estimate for pc- In particular 
one notes that the estimates obtained from /io,i(p) and /xo,2(p) (shown in the upper two 
panels on the right) do not appear to have converged yet. The estimates from the remaining 
series largely seem to have converged to a narrow band as the number of terms exceed 150. 
Secondly it appears that generally the scatter among the various estimates becomes smaller 
as the number of terms increases. All in all it appears that the estimates converge towards 
Pc — 0.644700185. The only notable exception is perhaps the estimates from the series 
l^2,o{j>)l^o,2ij>) I {.l^o,i{p)Y (shown in the lower left panel) which seems to favor a slightly larger 
value Pc — 0.644700191, though in this case there still seems to be a downwards trend in 
the estimates. It seems reasonable to estimate that pc = 0.644700185(5). This is in good 
agreement with the estimate Pc = 0.64470015(5) obtained previously |]2^. The new refined 
estimate is an order of magnitude more accurate and the central estimate lies within the 
error bounds of the earlier estimate. In Fig. ^ I have plotted estimates for the corresponding 
critical exponents. Similar trends are apparent in this case. From these plots I venture the 
following estimates 



7 


= 2.277730(5), 


^11 


= 1.733847(6), 


2u^ 


= 2.193708(4), 


7 + ^^11 


= 4.01156(1), 


7 + 2z/|| 


= 5.74539(1), 


7 + 2Z/X 


= 4.471425(15). 



(15) 



These estimates are again in full agreement with those obtained previously, an order of 
magnitude or so more accurate, and while the central estimates again are trending higher 
they still lie well within the earlier error bounds. Looking at differential approximants of 
different order (second and fourth) confirms the validity of these estimates. 

I also analysed the series in order to estimate the leading confluent exponents Ai. As 
was the case for the percolation probability series both the Baker-Hunter transformation 
and the method of Adler, Moshe and Privman (see [^ and references therein for details 



regarding these methods) yielded estimates consistent with Ai = 1. So there are no signs of 
non-analytic corrections to scaling. 

4.2 The site series 

An analysis similar to that of the previous section was performed for the site series. Generally 
I found that the estimates for pc were not so well converged. They seemed to change quite a 
bit as the order of the approximant or inhomogeneous polynomial increased, and furthermore 
the estimates also showed some inconsistencies from series to series. The series for fio,2{p) 
and /^o,i(p) yielded estimates for pc — 0.7054850 while the series /^2,o(p)/^o,2(p)/(/^o,i(p))^ 
favored an estimate Pc — 0.7054853 with the remaining series yielding estimates in between. 
From these one could surmise that Pc = 0.70548515(20). The associated estimate for the 
critical exponents are listed below 



7 = 2.27765(6), 
z/|| = 1.73381(4), 



2u^ 


= 2.19377(5), 


1 + n 


= 4.01135(15), 


1 + 2z/|| 


= 5.7450(2), 


7 + 2Z/X 


= 4.47130(8). 



(16) 



As can be seen these estimates are often only marginally consistent with those from the bond 
series and at least an order of magnitude less accurate. 

Due to the high degree of internal consistency of the estimates from the bond series 
one would tend to believe quite firmly in their accuracy and correctness and one can then 
use them to try and obtain a more accurate estimate for p^ in the site problem. In Fig. ^ 
I have plotted the estimates for the critical exponents vs the estimates for p^- The solid 
lines indicate the error-bounds on the exponent estimates obtained by using the bond-series 
estimates for the exponents 7, z/y and i/_l. By extrapolating the exponent estimates until they 
lie between the solid lines it is obvious that this procedure yields a pc estimate consistent 
with pc = 0.70548522(4), which I take as the final estimate. 

4.3 Non-physical singularities 

Non-physical singularities are of interest both because knowlegde about their position and 
associated exponents may help in the search for exact solutions and because one may gain a 
better understanding of the problem by studying the behaviour of various physical quantities 
as analytic functions of complex variables. While comparatively little work has been done 
along these lines for directed percolation this is a quite active field of study for classical 
spin systems such as the Potts and Ising models (see [3H| for recent results and references to 
earlier work) . 

The series have a radius of convergence smaller than pc due to singularities in the complex 
p-plane closer to the origin than the physical critical point. Since all the coefficients in 
the expansion are real, complex singularities always come in conjugate pairs. In order to 
locate the non-physical singularities in a systematic fashion I used the following procedure: 
Calculate all differential approximants with K and L fixed using at least 150 or 140 terms 
in the bond or site cases, respectively. Each approximant yields N^ possible singularities 
(and associated exponents) from the zeros of Qk (many of these are of course not actual 
singularities of the series). Next sort these 'singularities' into equivalence classes by the 
criterion that they lie at most a distance 2~'' apart. An equivalence class is accepted as a 
singularity if it contains more than Na approximants {Na can be adjusted but I typically 
used a value around 2/3 of the total number of approximants), and an estimate for the 
singularity and exponent is obtained by averaging over the approximants (the spread among 
the approximants is also calculated). This calculation was then repeated for fc — 1, k — 2, 
. . . until a minimal value of 6. To avoid outputting well converged singularities at every 
level, once an equivalence class has been accepted, the approximants which are members of 
it are removed, and the subsequent analysis is carried out only on the remaining data. This 
procedure is applied to each series in turn, producing tables of possible singularities. 

The analysis indicates that the series have quite a large number of non-physical singular- 
ities. A quick view of the distribution of singularities can be gained from Fig. |^ which show 
the location of the non-physical singularities. Estimates for the non-physical singularities 
are listed in Table |I|. For both the bond and site cases there are two sets of singularities. 
For those marked by 1, quite accurate estimates can be obtained for both the location of 
the singularity and the associated exponents. The exponent estimates for the site problem 
at the singularity on the negative axis are consistent with the exact values 1/2, —1/2, —3/2, 

in 



and —1/2 for the series S{p), /io,i, /io,2, and /i2,o, respectively. At the conjugate pair of 
complex singularities the corresponding exponents are consistent with the exact values 3, 2, 
1 and 2, respectively. In either case the error is no more than 0.1%. In the bond case the 
exponents at the conjugate pair of complex singularities seem to be identical to those for the 
site problem though the error is at least an order of magnitude larger. At the singularity 
on negative axis the exponent estimates are —0.075(5), —1.05(2), —2.025(15), and —1.00(2). 
So the most likely scenario for exact values would seem to be a logarithmic divergence in 
S(p), and divergences with exponents —1, —2, and —1 in /io,i, /io,2, and /i2,o, respectively. 
For the singularities marked 2, the estimates for the location of the singularity is poor and 
no meaningful estimates can be obtained for the exponents. In all cases the singularities are 
found in all series and remain fairly stable as K and L is varied. 

5 Conclusion 

In this paper I have reported on a new algorithm for the derivation of low-density series for 
moments of the pair-connectedness in directed percolation problems. Numerical evidence 
shows that the computational complexity grows exponentially, but with a growth factor 
A < v^, which is much smaller than the growth factor A = v^ of the previous best algorithm. 
For bond (site) percolation on the directed square lattice the series have been extended to 



order 171 (158) as compared to order 112 (106) obtained in previous work ^7_ . 

Analysis of the bond series led to a very accurate estimate for the critical point, pc = 
0.644700185(5), and the values of the critical exponents for the average cluster size, par- 
allel and perpendicular connectedness lengths are estimated to be 7 = 2.277730(5), z/y = 
1.733847(6) and z/_|_ = 1.096854(4). An accurate estimate for the percolation probability 
exponent is obtained from the scaling relation j3 = (z/y + z/^ — 7)/2 = 0.276486(8). 

For reference purposes I list estimates for some other critical exponents obtained using 
various scahng relations: 

A= /3 + 7= 2.554216(13) 

T= v\\- (3= 1.457362(14) 

z= z/||/z/^= 1.580745(10) 

7' = 7 - z^ii = 0.543883(11) 

5 = (3/v\\ = 0.159464(6) 

rj = 7/z/|| - 1 = 0.313686(8). 

Here A is the exponent characterizing the scale of the cluster size distribution, r is the cluster 
length exponent, z is the dynamical critical exponent, 7' the exponent characterizing the 
steady-state fluctuations of the order parameter, while 6 and r] characterize the behaviour 
at Pc as t — i> 00 of the survival probability and average number of particles, respectively. 

Assuming that the exponent estimates from the square bond case are correct, an improved 
Pc-estimate was obtained for the square site problem, pc = 0.70548522(4). 

Finally I note, that an analysis of the series, in order to determine the value of the 
confluent exponent, yielded estimates consistent with Ai ~ 1. Thus there is no evidence of 
non-analytic corrections to scaling. 
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E-mail or WWW retrieval of series 

The series for the directed percolation problems on the various lattices can be obtained via 
e-mail by sending a request to I.Jensen@ms.unimelb.edu.au or via the world wide web on 
the URL [ittp://www.ms.unimelb.edu.au/^iwan/| by following the instructions. 
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Table 1: Estimates of the location of various non-physical singularities in the series for 
directed bond and site percolation. 

Bond Site 

1 -0.51670(3) -0.451952165(8) 

1 -0.22605(10) ± 0.4400(l)i -0.2661783(4) ± 0.3847813(4)i 

2 0.135(10) ±0.455(10)i -0.0525(5) ± 0.4840(5)i 
2 0.0105(10) ±0.475(l)i 0.085(15) ± 0.500(15)i 
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Figure 1: The directed square lattice with orientation given by the arrows. 




Figure 2: Illustration of graphs which contribute to C{x,t) (solid hues), C{x',t') 
(solid+dashed lines) and paths which are not allowed (dotted lines). 
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Figure 3: Pictorial representation of the rules for updating boundary polynomials in the 
new algorithm. 
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Figure 4: The maximal number of contributing configurations in a calculation to order N. 
The straight line is a pure exponential oc A^ with growth factor A = \/2. 
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Figure 5: Estimates of the critical point Pc obtained from third-order differential approxi- 
mants vs the number of terms used by the approximant. Numbers along the y-axis are all 
preceded by 0.644700. Shown are (from left to right and top to bottom) estimates from the 
series S{p), /io,i(p), At2,o(p), Ato,2(p), At2,o(p)Ato,2(p)/(/io,l(p))^ and fJ-oM/l^oAp)- 
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Figure 6: Estimates of the critical exponents obtained from tliird-order differential approxi- 
mants vs the number of terms used by the approximant. Shown are (from left to right and top 
to bottom) estimates from the series S{p), /io,i(p), l^2,o{p), Aio,2(p), /^2,o(p)Aio,2(p)/(Aio,i(p))^, 
and/io,2(p)//io,i(p)- 
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Figure 7: Estimates of the critical exponents obtained from third-order differential approx- 
imants vs estimates of the critical point for square site problem. Numbers along the x-axis 
are all preceded by 0.70548. Shown are (from left to right and top to bottom) estimates 
from the series S{p), fio,i{p), At2,o(p), Ato,2(p), At2,o(p)Ato,2(p)/(Ato,l(p))^ and /io,2(p)//Uo,i(p)- 
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Figure 8: Location of non-physical singularities for bond (left panel) and site (right panel) 
directed percolation. 
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